Transmural APD heterogeneity determines ventricular arrhythmogenesis in LQT8 syndrome: Insights from Bidomain computational modeling

Long QT Syndrome type 8 (LQT8) is a cardiac arrhythmic disorder associated with Timothy Syndrome, stemming from mutations in the CACNA1C gene, particularly the G406R mutation. While prior studies hint at CACNA1C mutations’ role in ventricular arrhythmia genesis, the mechanisms, especially in G406R presence, are not fully understood. This computational study explores how the G406R mutation, causing increased transmural dispersion of repolarization, induces and sustains reentrant ventricular arrhythmias. Using three-dimensional numerical simulations on an idealized left-ventricular model, integrating the Bidomain equations with the ten Tusscher-Panfilov ionic model, we observe that G406R mutation with 11% and 50% heterozygosis significantly increases transmural dispersion of repolarization. During S1-S4 stimulation protocols, these gradients facilitate conduction blocks, triggering reentrant ventricular tachycardia. Sustained reentry pathways occur only with G406R mutation at 50% heterozygosis, while neglecting transmural heterogeneities of action potential duration prevents stable reentry, regardless of G406R mutation presence.

Previous experimental studies have shown that CACNA1C mutations can alter ICaL, potentially triggering ventricular arrhythmias.The recent in vivo study by Sanchez et al. [10] is particularly significant, as it investigates the arrhythmogenic mechanisms in the first knock-in swine model of TS1, demonstrating that cardiac activation impairments may lead to conduction blocks and reentry induction.
LQTS are linked to increased transmural and spatial dispersion of repolarization, a known precursor for reentrant ventricular tachycardia and fibrillation, [21,22].This dispersion is attributed to intrinsic differences in APD across the three principal cell types of the ventricular myocardium (endocardial, midmyocardial and epicardial), and to the extent to which these repolarization differences are damped by electrotonic forces, see [23][24][25].
Our study aims to elucidate the mechanisms by which the G406R CACNA1C mutation increases transmural dispersion of repolarization, potentially inducing and sustaining reentrant ventricular arrhythmia.We employ numerical simulations using the Bidomain model of electrocardiology [26][27][28][29], a non-linear reaction-diffusion system of Partial Differential Equations (PDEs), coupled with the ten Tusscher-Panfilov membrane model (TP06, see [30,31]), a stiff system of Ordinary Differential Equations (ODEs), describing the ion currents through the human ventricular cellular membrane.For numerical approximation, we use finite elements in space and semi-implicit finite differences in time, accelerated using parallel codes based on advanced Domain Decomposition techniques and run on Linux clusters with up to 128 cores.

Mathematical models
In the numerical section, we have simulated the bioelectrical activity of a three-dimensional idealized left ventricular geometry.We model the bioelectrical activity of a three-dimensional idealized left ventricular geometry using the Bidomain representation of the cardiac tissue [26][27][28][29].
Let O denote the three-dimensional region occupied by the cardiac tissue.Then, according to the Bidomain model, the evolution of the transmembrane potential v(x, t), extracellular potential u e , gating variables w(x, t) and ionic concentrations c(x, t) is described by the following system of partial differential equations: with appropriate initial conditions on v(x, 0), w(x, 0) and c(x, 0).Here c m and i ion denote the capacitance and the ionic current of the membrane per unit volume, i app represents the applied current per unit volume, D i,e are the intra-and extracellular transversely isotropic conductivity tensors.
Assuming transversely isotropic properties of the intra-and extracellular media, the conductivity tensors are given by where s i;e l ; s i;e t are the conductivity coefficients of the intra-and extracellular media measured along the fiber direction a l and any cross fiber direction, respectively.
The ionic current i ion (v, w, c) and the associated ordinary differential equations for the gating and ionic concentration variables are modeled according to the ten Tusscher-Panfilov membrane model (TP06) [30,31], available from the cellML depository (models.cellml.org/cellml), consisting of an ordinary differential equations (ODEs) system with 17 ODEs modeling the main ionic currents dynamics.The TP06 is one of the most used human ionic models with biophysical details and it has been extensively validated against experimental data from various experimental conditions, including different pacing rates, drug effects, and pathological states.Of course other biophysically detailed human ionic models could be considered as well, as long as they include the L-type calcium current (I CaL ) which in the LQT8 syndrome is altered by the CACNA1C mutation.

Numerical methods
The space discretization of the Bidomain system is performed by employing hexahedral isoparametric Q 1 finite elements, while the time discretization is based on the following double operator splitting procedure: a) split the ODEs from the PDEs and b) split the elliptic PDE from the parabolic one.The adopted numerical procedure has been described in detail in [32] and in the recent work [33].We also refer to [34][35][36] for other numerical strategies adopted in the literature.
For sake of clarity, we describe here the time discretization applied to the strong formulation of the Bidomain system.Let us denote by τ the time step size.At each time step, given v n , w n , c n , we proceed as follows: 1. find w n+1 , c n+1 by solving 2. find u n e by solving the elliptic PDE 3. find v n+1 by solving the parabolic PDE This operator splitting strategy yields at each time step the solution of two large scale linear systems of algebraic equations, arising from the finite element discretization of Eqs (3) and (4), respectively.
In order to ensure parallelization and portability of our Fortran code, we use the PETSc parallel library [37] developed at the Argonne National Laboratory.The two large linear systems at each time step are solved by a parallel conjugate gradient method, preconditioned by the Multilevel Additive Schwarz preconditioner, developed in [38], for the ill-conditioned elliptic system and the Block Jacobi preconditioner for the well conditioned parabolic system.These preconditioners are based on the multilevel PETSc objects PCMG (MultiGrid) with ILU(0) local solvers.The simulations are run on 64 cores of the Linux Cluster INDACO of the University of Milan (sito web) or on 128 cores the Linux cluster GALILEO of the CINECA laboratory (sito web).These computations are quite demanding: each run simulating 2000 ms after delivering a stimulus required about 12 hours on 64 cores, and for each of the 9 settings considered, we carried out several tests with 8 pacing runs and several S1-S4 stimulations as described in the Stimulation protocol section below.

Computational domain
The domain H is the image of a Cartesian periodic slab using ellipsoidal coordinates, yielding a truncated ellipsoid modeling a left ventricle (LV) geometry, described by the parametric equations where a(r , and a 1 = b 1 = 1.5, a 2 = b 2 = 2.7, c 1 = 4.4, c 2 = 5 (all in cm) and ϕ min = −π/2, ϕ max = 3π/2, θ min = −3π/8, θ max = π/8.We will refer to the inner surface of the truncated ellipsoid (r = 0) as endocardium and to the outer surface (r = 1) as epicardium.In all computations, a structured grid of 512 × 256 × 48 hexahedral isoparametric Q 1 finite elements of size h � 0.02 cm is used in space, for a total amount of 6447616 mesh nodes The time step size is fixed to τ = 0.05 ms.Fibers rotate transmurally, linearly with the depth and counterclockwise from epicardium to endocardium, for a total amount of 120 o .The values of the transversely isotropic conductivity coefficients are all given in mS.
In the TP06 model, the analytic expression of I * CaL , with � = WT, TS is where g CaL is the maximal conductance, Ca ss and Ca o the subspace and extracellular calcium concentrations, respectively, d* the activation gate, f* f 2 the voltage-dependent inactivation gates, f Ca ss the Ca ss -dependent inactivation gate, V the transmembrane potential, F the Faraday constant, R the gas constant, T the absolute temperature.All parameters values are those reported in [31].
On the basis of that, we model the ventricular tissue with homogeneous (HOM) or transmural heterogeneous distributions of APD.The transmural wall is subdivided into three layers of the same depth: the sub-endocardial layer (ENDO), the midmyocardial layer (MID) and the sub-epicardial layer (EPI).Two types of transmural heterogeneity are considered by scaling the original formulation (given in [31]) of the I Ks current by different factors: • in the HET1 setting, the I Ks current of the ENDO, MID and EPI cells, is multiplied by a factor 1.3, 0.5 and 1.4, respectively.As a result, in this heterogeneous setting, the midmyocardial cells present the longest APD, in agreement with the experimental studies [25,42,[45][46][47]; • in the HET2 setting, the I Ks current of the ENDO, MID and EPI cells is multiplied by a factor 0.7, 1.1 and 1.4, respectively.As a result, in this heterogeneous setting, the sub-endocardial cells present the longest APD, in agreement with the experimental studies [25,43,44].
In homogeneous simulation, all cells in the ventricle are assigned the calibration of ENDO cells, with I Ks multiplied by the factor 1.3.Fig 1 displays the resulting action potential waveforms, computed by solving the one-dimensional cable equation with conductivity coefficient 2 mS, whereas Table 1 reports the corresponding APD values at BCL = 1000 ms.Note that in case of TS-50%-HET1 (MID cells) and TS-50%-HET2 (ENDO cells), an EAD occurs.

Stimulation protocol
Our simulations use the S1-S2-S3-S4 stimulation protocol, a Programmed Ventricular Stimulation (PVS) which in synthesis is the pacing procedure introduced in Brugada et al. in 2003 [48] and frequently used in several experimental studies (see [49] for a recent review and alternative protocols in mice models).In this protocol, a first pacing cycle of several stimulations (S1) at a given cycle length is applied, and then three consecutive extrastimuli (S2, S3, S4) are applied at the closest coupling intervals at which each extrastimulus triggers an arrhythmia.
In our case, the first pacing cycle consists of S1 stimuli of 350mA/cm 3 for 1 ms on a small area 0.12×0.12×0.06cm 3 at an endocardial site for the HOM and HET1 setting and an epicardial site for the HET2 setting; the stimulus is delivered at a central location, equidistant from the apex and base of the idealized ventricle, as indicated by the breakthrough (red color) in Fig The S3 stimulus is initially delivered 350 ms after the previous S2 stimulus, and then shortened until arrhythmia is induced or the stimulus fails to induce arrhythmia.Analogously for S4.The criterion adopted for successful arrhythmia induction is the onset of a reentrant excitation that remains sustained up to the end of simulation time, set to 4000 ms after the last stimulus.This stimulation protocol is quite demanding computationally, since each run simulating 2000 ms after delivering a stimulus requires about 12 hours on 64 cores of our Linux clusters.This cost, scaled by the simulated time, occurs for each of the 8 pacing runs and S1-S2-S3-S4 stimulations, which we repeated in each of the 9 settings considered, from WT-HOM to TS-50%-HET2 (see previous section).

Post-processing
In each simulation, we save the time evolution of the transmembrane potential on a matrix of 64 × 33 exploring multi-electrode transmural needle, equally distributed through the LV.Each needle carries 7 recording sites, equally spaced along the transmural epi-to-endocardial direction.Therefore, we have a total of 14784 recording sites.
From each transmembrane potential waveform, we compute: • the activation time (ACTI), defined as the instant of maximal slope during the depolarization phase of the action potential; • the repolarization time (REPO), defined as the instant of minimal slope during the repolarization phase of the action potential; • the action potential duration (APD), defined as the difference APD = REPO-ACTI.
As a byproduct, we compute the following quantities: • total activation time dispersion (dAT tot ), defined as dAT tot = max(ACTI)-min(ACTI), where max and min are taken over all the 14784 recording sites; • total repolarization time dispersion (dRT tot ), defined as dRT tot = max(REPO)-min(REPO), where max and min are taken over all the 14784 recording sites; • average transmural repolarization dispersion (dRT trans ), defined as the average over all the 64 × 33 transmural needles of max(REPO)-min(REPO), where max and min are taken over the 7 sites of each needle; • total APD dispersion (dAPD tot ), defined as dAPD tot = max(APD)-min(APD), where max and min are taken over all the 14784 recording sites; • average transmural APD dispersion (dAPD trans ), defined as the average over all the 64 × 33 transmural needles of max(APD)-min(APD), where max and min are taken over the 7 sites of each needle

Numerical results
We have simulated the electrical activity of the idealized LV model, considering the six settings described above in the Simulation setup section.For each setting, we first perform 8 beats at a BCL = 500 ms.Then we apply an S1-S2-S3-S4 stimulation protocol by reducing the coupling interval between subsequent stimulations.In HOM and HET1 simulations, the stimulus is applied at an endocardial site, whereas in HET2 simulations the stimulus is applied at an epicardial site.

Total and transmural dispersion of ACTI, REPO, APD
In Table 2, we report the dAT tot , dRT tot , dRT trans , dAPD tot and dAPD trans parameters, as computed after the S1 stimulus, thus corresponding to a BCL = 500 ms.In the settings without transmural heterogeneities (WT-HOM, TS-11%-HOM, TS-50%-HOM), we do not observe significant differences between the pathological settings (TS-11% and TS-50%) and the WT setting, in terms of all five parameters.On the other hand, in the settings with transmural heterogeneities, we do observe a strong increase of total and transmural dispersion of repolarization and APD in case of TS-50%-HET1 with respect to WT-HET1 (331 vs 195 ms for dRT tot , 94 vs 25 ms for dRT trans , 144 vs 39 ms for dAPD tot , 97 vs 27 ms for dAPD trans ) and in case of TS-50%-HET2 with respect to WT-HET2 (382 vs 218 ms for dRT tot , 145 vs 32 ms for dRT trans , 196 vs 42 ms for dAPD tot , 142 vs 31 ms for dAPD trans ).
The results have also shown that, in the settings without transmural heterogeneities, the average APD on the entire LV increases from 267 ms (WT-HOM) to 286 ms (TS-11%-HOM) and 346 ms (TS-50%-HOM), whereas, in the settings with transmural heterogeneities, it increases from 277 ms (WT-HET) to 297 ms (TS-11%-HET) and 408 ms (TS-50%-HET).This latter increase is mainly due to the increase of APD in the midmyocardial layers.
In order to better understand the dynamics of the activation and repolarization processes, Figs 2 and 3, display the spatial distributions of ACTI, REPO, APD, computed after the S1 stimulus, on the epicardial, midmyocardial, endocardial and on two longitudinal transmural sections, for the WT-HET1 and TS-50%-HET1 settings.Note that, on the midmyocardial section, TS-50%-HET1 presents a strong increase of REPO and APD dispersion with respect to WT-HET1 (240 vs 158 ms for REPO, 85 vs 17 ms for APD).These dispersions are computed as the difference between the maximum and minimum value of the corresponding distribution (REPO or APD), on the midmyocardial section.Analogously, Figs 4 and 5, report the spatial distributions of ACTI, REPO, APD for the WT-HET2 and TS-50%-HET2 settings.In this case, on the endocardial section, TS-50%-HET2 presents a strong increase of REPO and APD dispersion with respect to WT-HET2 (252 vs 158 ms for REPO, 189 vs 16 ms for APD).

Arrhythmogenesis with S1-S2-S3-S4 stimulation protocol
In the TS-11%-HET1 setting, after the S1 stimulus, we apply an S2 stimulus with an S1-S2 coupling interval of 320 ms, and an S3 stimulus with an S2-S3 coupling interval of 260 ms.After the S4 stimulation delivered at an S3-S4 coupling interval of 240 ms, a conduction block occurs when the excitation wavefront reaches the midmyocardial region, because this latter region presents a larger APD and thus effective refractory period than the endo-and epicardial regions.Then the electrical signal proceeds around the region of block, towards epicardium.After having excited the epicardial surface, the excitation wavefront comes back into the midmyocardial regions, now excitable again, generating a reentrant ventricular tachycardia, which dies after three cycles of reentry, as confirmed by the transmembrane potential and unipolar electrogram waveforms reported in Fig 6 (see also S1 Movie in the supplementary data).In the TS-50%-HET2 setting, after the S1 stimulus, we apply an S2 stimulus with an S1-S2 coupling interval of 360 ms.After the S3 stimulation delivered at an S2-S3 coupling interval of 340 ms (Fig 8(a)), a conduction block occurs when the excitation wavefront reaches the subendocardial region, which presents a larger APD and effective refractory period than the midmyocardial and sub-epicardial regions.Then the electrical signal proceeds around the region of block, towards the endocardium.After having excited the endocardial surface, the excitation wavefront comes back into the midmyocardial regions, now excitable again, yielding an epicardial BKT (Fig 8(c)), at about 360 ms after the S3 stimulus, and thus triggering the first cycle of reentrant excitation.After having excited the epicardium, the reentrant wavefront moves towards the endocardium, reactivating first the apical regions.The activation wavefront moves again back to epicardium, yielding two BKTs ( In the WT-HET1, WT-HET2, TS-11%-HET2 settings and in all the settings without transmural heterogeneities, no conduction blocks occur during the entire S1-S2-S3-S4 stimulation protocol.Consequently reentry is not induced in these cases.

Discussion
We have investigated by computer modeling the role played by transmural APD heterogeneity on the onset of reentrant ventricular arrhythmias in the presence of LQT8 syndrome.This congenital heart disease might be produced by several mutations in the CACNA1C gene, encoding the CaV1.2 subunit of the ICaL current.The mutation considered here is the G406R mutation [7,8].
The three-dimensional numerical simulations reported are based on the Bidomain model of electrocardiology, coupled with the TP06 membrane model, taking into account the fiber anisotropy of the ventricular tissue.Our calibration of the G406R mutation in the ICaL current provided by the TP06 model is in agreement with the experimental and computational findings reported in [8].
We have considered two types of transmural APD heterogeneity: in the first one (HET1), midmyocardial cells (M cells) present a prolonged APD with respect to sub-endocardial and sub-epicardial cells; in the second one (HET2), sub-endocardial cells present a prolonged APD with respect to midmyocardial and sub-epicardial cells.
The results have shown that, with both HET1 and HET2 APD heterogeneities, the presence of G406R mutation with 50% of heterozygosis (TS-50%) yields a significant increase of transmural dispersion of repolarization.This growth is due to the prolonged APD of the midmyocardial and sub-endocardial cells in the HET1 and HET2 settings, respectively.
During the implementation of the S1-S2-S3-S4 stimulation protocol, the high transmural repolarization gradients occurring in the TS-50%-HET1 and TS-50%-HET2 simulations cause conduction blocks of the excitation wavefronts elicited by premature beats.Such blocks of conduction induce the onset of reentrant ventricular tachycardia.EADs occuring in the regions of prolonged APD facilitate the maintenance of the cycles of reentry, as observed in previous computational studies [50,51] for other LQT variants.When transmural APD heterogeneities are neglected, in WT-HOM, TS-11%-HOM and TS-50%-HOM simulations, reentry does not occur, irrespective of the presence of LQT8 pathology.The onset of reentry due to dispersion of repolarization after premature beats and the occurrence of EADs was also observed in the recent experimental study [10].However, in [10], the mechanism responsible for the increase of dispersion of repolarization was an activation delay due to impairment of I Na current, rather than the combination of G406R mutation and intrinsic APD heterogeneities.
We recall that in our HET1 heterogeneous setting, midmyocardial cells (M cells) present a prolonged APD.In vitro experimental studies have demonstrated that isolated myocytes, extracted from different depths across the ventricular wall, exhibit transmembrane action potentials with different durations, see [40,41], with the APD of M cell being longer than that of sub-endocadial and sub-epicardial cells.This transmural APD heterogeneity has also been observed in myocardial wedge preparations, see [25,42,[45][46][47].Despite these findings obtained in myocardial in vitro and wedge preparations, high transmural repolarization gradients have never been observed in intact hearts, see e.g.[52][53][54][55].It has been suggested that the electrotonic coupling present in vivo mask cellular heterogeneities, yielding a reduced dispersion of repolarization times and APDs.Indeed in our WT-HET1 simulation, the resulting transmural dispersion of repolarization time and APD is not large (25-32 ms for repolarization time and 39-42 ms for APD), and it is comparable with those observed in intact hearts.
In our HET2 heterogeneous setting instead, an epi-to-endo transmural APD gradient is introduced, with sub-endocardial cells presenting a prolonged APD.This setting is in agreement with recent experimental studies on human wedge preparations reported in [25,43,44].

Limitations
In our model of I CaL current, we have assumed that the G406R mutation affects only the voltage dependent inactivation gate.However, experiments reported in the recent study [39] have shown that the G406R mutation might affect also the calcium dependent inactivation gate of the I CaL current.Future tissue simulation studies should investigate how the G406R mutation in the I CaL calcium dependent inactivation gate influences the onset of arrhythmic events.Moreover, a future computational study should also take into account the reduction of I Na current, as suggested in the experimental work [10].
for 50% channels open V * ina0:5 and scaling parameter S * ina given by • V WT ina0:5 ¼ À 20:0 mV, S WT ina ¼ 7:0; • V TS ina0:5 ¼ 0:0 mV, S TS ina ¼ 12:0.The parameters were adjusted so that our activation (d * 1 ) and inactivation (f * 1 ) curves, displayed in Fig 1, match at least qualitatively the experimental plots reported in Fig 4E and 4F of [8], and the resulting effects in our simulated pathological action potential waveforms, also displayed in Fig 1, are comparable with those observed in Fig 6 of [8].The waveforms in Fig 1 are computed by solving a one-dimensional version of system (1).

Table 1 . Action potential duration of WT and TS endo-, mid-and epicarical cells at BCL = 1000 ms.
7, panel a) at t = 380 ms, and in Fig 8, panel a) at t = 340 ms.For each simulation setting, we first apply 8 pacing stimuli (S1) at a BCL of 500 ms.Then, a premature stimulus (S2) is delivered 380 ms after S1.If S2 does not generate a reentrant arrhythmia, the S1-S2 coupling interval is shortened in steps of 10 ms until arrhythmia is induced or S2 fails to trigger excitation.If arrhythmia is not induced, an additional S3, and if necessary, an additional S4 stimulus, is delivered in the same manner as S2.